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DIFFUSION IN SPATIALLY VARYING POROUS MEDIA 

MARIA BRUNAtt AND S. JONATHAN CHAPMANt 


Abstract. The problem of diffusion in a porous medium with a spatially varying porosity is 
considered. The particular microstructure analyzed comprises a collection of impenetrable spheres, 
though the methods developed are general. Two different approaches for calculating the effective 
diffusion coefficient as a function of the microstructure are presented. The first is a deterministic 
approach based on the method of multiple scales; the second is a stochastic approach for small volume 
fraction of spheres based on matched asymptotic expansions. We compare the two approaches, and 
we show good agreement between them in a number of example configurations. 
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1. Introduction. The macroscopic transport of solutes in porous media de¬ 
pends critically on microscopic features such as the structure of the porous matrix 
and the nature of the interactions between the solid and liquid phases. On the other 
hand, the complexity of the microscopic problem means that in practice it is often 
desirable to obtain a macroscopic effective-medium equation from which the macro¬ 
scopic transport can be obtained directly. This idea of upscaling is ubiquitous in 
many sciences. In particular, diffusive transport in heterogeneous media occurs in hy¬ 
drogeology (aquifers, groundwater [13, 14]), contaminant transport (water filtration 
with membranes), lithium-ion batteries [7], and biological applications such as porous 
biofilms [12] and intracellular transport [22]. 

Efforts to find the form of the macroscopic equation and estimate the effective 
properties of the medium date back at least to the 19th century, when Maxwell gives 
an approximate expression for the effective conductivity of a heterogeneous medium 
comprising small (i.e. well-separated) spheres of one material distributed in another 
[21, p.403]. Since then a great variety of approaches have been developed to obtain 
such upscaled equations [9]. A usual starting point is to suppose that the solute 
undergoes a simple diffusion process in the void or fluid phase of the porous 
medium 11 so that the evolution of the concentration of particles C'(x, t) is described 

by 

dC 

— =V-(UoVC), xeO„, (l.la) 

ot 

0 = n-(UoVC'), xeail„, (1.1b) 

where d^y is the solid-fluid interface and n is the outward unit normal to 11„. The 
fraction of space available to the diffusing species is the porosity dr = |n„|/|H|. If 
the solid phase is denoted by Us (the complement of n„) then the volume fraction of 
solid is $ = 1 — Ik. When the microstructure is fine we would like to be able to use 
an effective transport equation such as 

3c 

— = V-[UeVc], xell = ll,ull„ (1.2) 
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where -Dg is an effective diffusion coefficient, which will depend both on Dq and the 
geometry of Here c is the homogenized solute concentration (whose definition will 
be made precise in the next section). Note that c is defined throughout the material, 
whereas C is defined only in the fluid region. 

In general, the homogenization procedure will depend on both the porous medium 
structure and the transport processes within the medium. A large number of homog¬ 
enization techniques have been developed over the years [9]. These can be divided 
into two very broad categories, deterministic and stochastic techniques. Deterministic 
techniques include techniques such as volume averaging [34], multiple scales [29] and 
the generalized Taylor-Aris-Brenner method of moments [3]. Here the averaging relies 
on the separation of scales between the microstructure and the macroscopic material, 
and is a local spatial average. Many techniques, included the method of multiple 
scales and the generalized Taylor-Aris-Brenner method of moments, assume some 
periodicity in the microstructure. In this case, the microscopic lengthscale measures 
variations within a period cell and the macroscopic lengthscale measures variations 
within the macroscopic region of interest. More specifically, the method of multiple 
scales constructs a family of problems involving the ratio 5 between these two length- 
scales and uses asymptotic expansions in 5 to study systematically the convergence 
as (5 —?► 0 to a limit problem. Since 5 is assumed to be small, the resulting family of 
problems contains rapidly oscillating (periodic) coefficients. Thanks to its systematic 
nature, this method has formed the basis of a whole field in mathematics known as 
mathematical homogenisation [1]. 

In terms of the example above, deterministic approaches upscale from (1.1) to 
(1.2) by assuming a particular given geometry of the solid matrix VLs from which the 
effective coefficient De can be computed. For example, it is common to consider either 
very simple periodic structures [such as the array of spheres depicted in Figure 1.1(a)] 
or to consider a disordered unit cell representative of the material as a whole, which is 
then extended periodically [31]. In the case of a periodic array of spherical inclusions, 
Df, can be computed exactly using Rayleigh’s multipole method [26], leading to an 
infinite system of algebraic equations. More complex microscopic structures can be 
dealt with via multiple scales (see for example [27]) or volume averaging [34]. Although 
these two methods are based on different underlying principles, they result in the 
same averaged equations [11]. These equations contain effective parameters (such 
as De) which depend on the microscale and are evaluated using the solution of a 
cell problem in the periodic unit cell (in multiple scales) or the solution of a closure 
problem in a representative elementary volume (in volume averaging). For a review of 
the application of the method of volume averaging in ordered and disordered porous 
media we refer the reader to [25]. 

The assumption of periodicity can be seen as artificial or too idealistic when 
modelling real heterogeneous media. As a result, an obvious progression came in the 
form of homogenization techniques for random media, in an attempt to reflect the 
uncertainty caused by the high degree of heterogeneity as well as the lack of experi¬ 
mental data [10]. In its simplest setting, this can be seen as replacing the periodicity 
assumption in the multiple scales method by stochastic periodicity (or statistical ho¬ 
mogeneity) [19, 23]. In particular, this means that the problem now contains random 
coefficients rather than periodic ones, and the homogenization consists of taking an 
averaging window of size <5 large enough so that ergodicity holds, that is, that a spa¬ 
tial average is equivalent to an ensemble average. A rigorous derivation of the limit 
problem as d —>■ 0 is quite challenging in general, but has been done in specific cases 
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such as steady heat conduction [19, 23] or, in a discrete setting, in a random con¬ 
ductance model for a network of resistors [16]. In natural porous media such as soils 
and aquifers, the problem to be upscaled is usually more challenging than the cases 
above, as often one must consider the coupled problems for groundwater flow and so¬ 
lute transport. Stochastic approaches developed in the context of hydrology include 
the works of Cushman [8], Dagan [10], and Gelhar [15]. A common starting point 
is to take the hydraulic conductivity to be a random field, for example assuming it 
has a lognormal probability density function [10] (p. 164). In particular, this implies 
that the Stokes velocity and the solute concentration are random fields as well. The 
classical approach is to linearize both the flow and the transport equations and to 
then take ensemble averages of the resulting problem [9]. 

In terms of the particular problem (l.I), a natural stochastic approach is to 
suppose that the microstructure is random and statistically homogeneous, with the 
effective properties arising by taking an ensemble average over different realizations 
of the microstructure [31]. For example, one could assume that the solid matrix fig 
is composed of spherical inclusions uniformly distributed in fl with a non-overlapping 
constraint [see Figure 1.1(b) for one sample of such configuration]. However, it is 
usually very difficult to estimate the effective properties from such a description, 
since doing so requires an infinite set of statistical correlation functions, which are 
generally never known [31]. Brown [4] provided a series form of the effective constant 
De (for dielectric constants rather than diffusivities, but the problem is analogous) 
and showed that the first correction term to Dq depends only on the solid volume 
fraction $, but that all other terms require knowledge of the statistical distribution 
of the solid matrix. Thus the best estimate one can obtain knowing only $ and Dq is 
equivalent to Maxwell’s formula, namely De = Do/{l + ^/2)* As a result, subsequent 
work by Hashin and Shtrikman [17] focused on using variational methods to obtain 
lower and upper bounds on De as a function of $, independent of the statistics of 
the medium. Prager [24] developed a method to introduce two- and three-particle 
correlation functions into such bounds. A specific application of these variational 
principles in the case of a porous medium formed by uniformly distributed spheres 
without non-overlapping constraints was used by Weissberg [33], who obtained the 
upper bound De < Do/{l — ^ ln(l — $)).! 

It is not completely clear which of these approaches is preferable in any given sit¬ 
uation, or how they compare to each other. As pointed out in [20], enforced periodic 
structures can have an important influence on the homogenized model in low-porosity 
systems. On the other hand, the variational approaches used by Brown [4] or Weiss¬ 
berg [33] are very general and do not require any constraints on scale or periodicity. 
They would be applicable also for so-called non-homogenizable media, in which there 
is no clear separation of scales between micro and macrostructures. 

One situation in which all the approaches discussed above need to be rethought is 
that in which there is macroscopic inhomogeneity in the material as well as microscopic 
inhomogeneity. Such a situation would arise in a material in which there is a non- 
uniform porosity on the macroscale. It might be thought that since the porosity is 
varying slowly it is locally constant, and the standard formulae can be used albeit 


*The effective diffusivity reported was in fact De = De{l — '!>) in the current notation. Thus, the 
expression reported in [33] for example is De/Do = (1 — + ‘I’/2). The two definitions result in 

the same macroscopic equation for homogeneous media, but not when the volume fraction varies on 
the macroscale, as we will see later. 

t Again, this expression is reported as a volume averaged diffusivity, De = De(l — <I>) in [33]. 
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Fig. 1.1. Two examples of porous media, (a) Periodic structure (with spherical inclusions 
on a square lattice), (b) One sample from a random structure, with spherical inclusions uniformly 
distributed and with non-overlapping constraints. Both media correspond to a solid volume fraction 
0 /$ = 20 %. 


with a cell problem (and hence diffusivity) which varies on the slow scale. However, 
we will see that this is not the case. 

Despite its acknowledged importance in many physical systems (for example, in 
the design of membranes for water filtration, in groundwater flow systems [13], and in 
Portland cement [28]), non-uniform porosity has rarely been accounted for in models 
and the effect on diffusion is not well understood. One exception is [32], which uses 
the method of volume averaging in systems with gradients in porosity. Worryingly, 
though, they find that the effective diffusivity, the porosity, and its gradient are highly 
dependent upon the location of the centroid of the representative elementary volume. 

In this paper we derive the effective transport equation for diffusion in a porous 
medium with gradients in porosity. We consider two types of porous media, namely an 
ordered medium with quasi-periodic inclusions and a disordered medium composed of 
randomly (but not uniformly) distributed inclusions. The homogenization procedures 
for each of these media are fundamentally different. Thus, the purpose of the paper 
is twofold: first, to study how macroscopic changes in the microstructure affect the 
homogenized equation and, second, to compare the two different approaches used to 
model diffusion in porous media. 

We consider a simple diffusion process in the void or liquid phase, where the solute 
has a constant diffusion coefficient Dq. The solid phase is composed of spherical and 
non-overlapping inclusions, which are impenetrable by the solute (that is, the diffusion 
coefficient is zero in the solid phase). For most of this paper we consider a fixed 
solid matrix (although we later remove this assumption in the stochastic approach to 
consider moving obstructions). 

We obtain the effective diffusion coefficient as a function of the microstructure of 
the medium, and in particular the porosity. We find that a non-uniform porosity on 
the macroscale results in an advection term in the homogenized equation, as found 
by [32]. This term accounts for the biasing of the motion of solute particles towards 
regions of high porosity. We also show that our two approaches, although very dif¬ 
ferent, can be reconciled to give the same macroscopic equation in the limit of small 
solid volume fraction. 

The article is organized as follows. In §2 we consider a deterministic approach 
and we apply the method of multiple scales to derive the effective transport equation 
in ordered porous media. We show how the method can be adapted to non-uniform 
porosity and non-periodic structures, providing they are locally periodic. We then 
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consider the same problem for stochastic porous media (disordered and non-periodic) 
in §3, and derive an equivalent effective equation from the Fokker-Planck description 
of the microstructure. In §4 we present two numerical studies coirrparing the two 
approaches, both with each other and with the exact solution of the corresponding 
microscopic problems. Finally, in §5, we present our conclusions. 

2. An ordered porous medium. We rescale length and time so that the 
rescaled domain has volume one, |f2| = 1 and the molecular diffusion coefficient is 
Do = 1.^ We also rescale concentration so that Cdx = 1. For simplicity we as¬ 
sume that fi = [—1/2,1/2]'^, where the dimension d may be 2 or 3. We suppose that 
the solid phase consists of Ng spherical and non-overlapping obstacles with radii and 
centers r^ for i = 1,... Ng. We suppose that the radius is a slowly varying function 
of position given by e(x), so that we may write = e(ri). Thus ftg = U^\i3£(r.)(ri), 
where Bg{r) is the d-dimensional ball of radius e centered at r. 

We begin by supposing that the centers r^ lie on a regular square lattice with 
period d <C 1, as depicted in Figure 2.1. Then (1.1) reduces to 

dC 

— = V^C, X G n,, (2.1a) 

u-VC' = 0, on ||x — Till = e(ri), l<i<Ng, (2.1b) 

with Ng = S~‘^. The local porosity is then iIj{x) = 1 — The global 

porosity is the average of the local porosities and is determined by the volume of all 
Ng obstructions, namely 


. N, Ns 

® i=l i=l 

2.1. Asymptotic homogeuizatiou via the method of multiple scales. We 

use the method of multiple scales to derive an averaged (or homogenized) model for 
the concentration C, valid over many obstacles, in the limit d <g; 1. We retain x as 
the macroscale variable, measuring distance on the scale of the whole sample, and 
introduce y = x/d as the microscale variable, which measures distance over the scale 
of the obstacle separation. As is usual in the method of multiple scales we treat these 
two variables as independent. The extra freedom this gives is removed by enforcing 
that the solution is exactly periodic in y; small variations from one unit cell to the 
next are thereby captured through the macroscale variable x. 

After introducing these two scales and using the chain rule spatial derivatives in 
(2.1) transform according to 


V-AVx + ^Vy. (2.2) 

We introduce e(x) = e{x)/S as the the obstacle radius relative to the dimension of the 
unit cell. As stated above, a crucial assumption we make is that this varies slowly, that 
is, it depends on x but not on y. We denote the unit cell by y G uj = [—1/2,1/2] 
and the solid phase by y G a;s(x) = i?£(x)(0)i with the fluid phase then given by 
y G w„(x) = uj\ 0Jg{x) (see right-hand side of Figure 2.1). 

^We abuse the notation by using the same symbols D and Do as in the dimensional problem 
( 1 . 1 ). 
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X G f2 y C w 




Fig. 2.1. Illustration of the problem geometry in two dimensions. Left: macroscale domain 12, 
with obstacles of radius e(x) placed in a periodic array. Right: microscale domain u), unit cell with 
an obstacle of radius e(x) = e(x)/5. 


The expression (2.2) allows us to transform (2.1a) according to the multiple-scales 
approximation. What is not clear is how we should write (2.1b) in multiple-scales 
form, given that the radius of the obstacle depends on x, so that the unit normal 
depends on x as well as y. At first sight it would seem that geometry dictates that 


n = ny = 


y 

l|y|| 


(2.3) 


and (2.1b) should be written 

d- Vx^ C • Hy = 0 for y G dujy{x). (2.4) 

However, this is incorrect as it neglects the variation of n with x. 

The systematic way to derive the multiple-scales equivalent of (2.1b) is to intro¬ 
duce the function 


X(x,y) = e(x) - ||y||, (2.5) 

for which the fluid-solid interface is the level set x(x, y) = 0. Note that this idea can 
be extended to more complicated interfaces in a straightforward manner. The normal 
is in the direction Vy, which we can now readily put in multiple-scales form using 
( 2 . 2 ): 

noc VxX + ^VyX = Vxe(x)-i||^. (2.6) 

This expression for the normal in multiple-scales form is far from obvious, and is a 
crucial part in extending the theory in [27] to macroscopically varying cell geometries. 
There are two contributions to the interface normal in (2.6); the first is due to the 
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slow variation of the obstacle radius, and the second is the unit radial vector in the 
microscale. 

Thus writing (2.1) in multiple-scales form gives 


dC 

(5^— = Vie + (5Vx • VyC -b (5Vy • VxC -b 6‘^VlC 

dt y y y 

VyC ■ ny = -SV^C ■ ny - SVyC ■ Vx£(x) - S^V^C ■ Vxe(x) 


yGa;.„(x), (2.7a) 

y € i9w„(x), (2.7b) 


where ny (given in (2.3)) is the unit normal vector into the obstacle in the microscale 
variable. 

Our aim is to derive the equation for the averaged macroscopic concentration 
c(x, t) over a representative volume located at position x. With our periodic mi¬ 
crostructure we may define this as 


c(x,t) 


tK [ C'(x,y,<)dy 
1^1 Joj 




C'(x,y,t)dy, 


( 2 . 8 ) 


since C = 0 in the solid phase. In the context of volume averaging, the unit cell ui is 
the representative elementary volume (REV), and is located at the centre of each of 
the obstacles. The average c is referred to as volumetric or superficial average. It is 
related to the intrinsic average 


c(x,t) = -— C(x,y,t)dy, 


through the porosity '!/’(^)- 


c(x, t) ='(/'(x)c(x, t) where 'tp{x) = 


;(x)| 


w 


(2.9) 


( 2 . 10 ) 


With a spherical obstruction in the centre of each cell, the porosity is = 1 ~ 

|i3g(x)|. The concentration c is normalized in the whole space, whereas c is normalized 
in the available space. 

Following the standard multiple-scales method, we now seek a solution to (2.7) in 
the limit of small <5 of the form C = C{x,y,t) which is periodic in y, while treating 
X and y as independent. Expanding in powers of 6 as C{x,y,t) = C(°^(x, y, t) -b 
6C^^\x,y,t) -b S^C^‘^\x,y,t) + ■■■ we find the leading-order equations require that 
C(o) is independent of y. At first order in S we find 


= 0 

y G a;„(x). 

(2.11a) 

VyC(^) • ny = -VxC(°) • ny 

y G dujy{x), 

(2.11b) 

periodic 

in y. 

(2.11c) 

The solution of (2.11) can be written as 



CW(x,y,t) =-VxC(°)(x, 

t) • r(x,y). 

(2.12) 


where r(x,y) is a d-vectorial function, whose components T^ satisfy the following cell 
problem: 


= 0 

VyTj ■ Ily = Tiyi 

Fi periodic 


y G w^(x), 
y G dujyix), 

in y, 


(2.13a) 

(2.13b) 

(2.13c) 
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where Uy^i is the ith component of the unit vector ny. We note that F varies with 
the macroscale variable x because of the variation of with x. 

Proceeding to second order in the expansion of (2.7) leads to the following problem 
for (7(2): 

= Vx • (VyCd) + VxC'C’)) + Vy • (VyC(2) + VxC(i)) Y S (x), 

(2.14a) 

VyC'(2) • ny = -Vx(7(i( • ny - (Vy(7(i( + Vx(7(°)) • Vxe(x) y G (9w„(x), 

(2.14b) 

( 7 ( 2 ) periodic in y. (2.14c) 

Integrating (2.14a) over a;„(x) using the divergence theorem, (2.14b) and periodic 
boundary conditions on (7(^) and (7(2) yields 


'0(x) 


9(7(0) 

dt 




(Vy(7(i( + Vx(7(°() dy 

Vxe(x) d5'y. 


[ (VyC(l)+V: 


(2.15) 


Now using the transport theorem to switch the order of integration with respect to y 
and differentiation with respect to x in the term on the right-hand side gives 


V'(x) 


a(7(°) 

dt 



(2.16) 


Let Jr(x,y) be the Jacobian matrix of F given by (Jr)ij = dVi/dyj. Then using 
(2.12) we may write 

VyC(i) = - (x, y) VxC'(°(. (2.17) 


Combining (2.16) and (2.17) we arrive at the following equation for the leading-order 
intrinsic average c = (7(°) (x, t): 

dc 

•0(x)— = Vx • ['0(x)Z)e(x)Vxc], (2.18a) 

where the effective diffusion tensor is given by 

He(x) =/d - —^ / Jf(x,y)dy, (2.18b) 

V(x) ^4x) 

where Id is the identity matrix of dimension d. The resulting model (2.18) is equivalent 
to that derived by [32] using a volume averaging approach. § 

We see that if ijj is independent of x then it may be cancelled on both sides of 
(2.18a) giving the usual homogenized result from multiple scales. However, when the 
radius of the obstacles varies with x the effect is felt not only through De (due to a 
cell problem that depends on x) but also through the fact that ip on the right-hand 
side appears inside the x-derivative. 


§See their equations (16) and (17); in our notation, = blx), (Cyi-y)"’' = c{'x.,t) and Deff = 
Z)e(x). However, we note that their cell problem (14) includes formally higher-order terms and 
hence is not exactly equal to our cell problem (2.13). 
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2.1.1. Comparison between averages: superficial versus intrinsic. The 

homogenized model (2.18) for the intrinsic average c is not yet an effective diffusion 
equation because of the porosity ip multiplying the left-hand side of (2.18a). Hence 
it is incorrect to refer to '0(x)lle(x) in (2.18a) as an effective diffusion coefficient. 
However, if we rewrite it in terms of the volume average c = tjjc we find 


m 


= Vx- 





= Vx- 


£'e(x)VxC - 


£>e(x)Vx'0(x) ■ 

'0(x) 


(2.19) 


Thus, we see that, when written in terms of c, the equation turns into a classical 
advection-diffusion equation, with Ile(x) the effective diffusion coefficient. The ad- 
vection term v = Z)e(x)Vx'0(x)/^/^(x) arises when the porosity is non-uniform. We 
see that solutes diffuse with a bias towards regions of increased porosity. 

In the case of uniform porosity '0(x) = dt, the intrinsic average c is simply a 
scaled version of the volumetric average c and their respective equations (2.19) and 
(2.18a) both reduce to a diffusion equation with constant diffusion coefficient as 
in (1.2). However, when the porosity is non-uniform the equations are different, and 
it is important to be clear which average of C we are dealing with. 

2.1.2. The diffusion tensor. To evaluate the diffusion tensor Df,{x) in (2.18b), 
we must solve the cell problem (2.13). So far we have not used that the obstacle in 
each cell is spherical, and in principle we could solve the cell problem (numerically) 
for any obstacle shape. However, having a spherical obstacle greatly simplifies the 
procedure. This is because, by symmetry, we find that is a multiple of the identity 
tensor, so that we have a single scalar diffusion coefficient. 

We solve the cell problem (2.13) using COMSOL Multiphysics and evaluate the 
integrals in (2.18b) numerically. We repeat for various 0 < e < 0.5 and plot the 
resulting Dg in Figure 2.2(a) as a function of the solid volume fraction (/)(x) = 7re'^(x). 
The three-dimensional counterpart is shown in Figure 2.2(b). 




Fig. 2.2. Effective diffusion versus volume fraction (j>. Shown are the simulation results of Dg 
in (2.18b) (data points), Rayleigh’s multipole method values (solid lines) and asymptotic result De 
(3.7a) (dash lines), (a) d = 2; Square lattice [square data points and (2.20a)/ and hexagonal lattice 
(circle data points and (2.20b )); asymptotic value Dg ~ 1 — 0. (b) d = 3: Cubic lattice (square data 
points and (2.20c)/ and asymptotic value De ~ 1 — <//2. 

We next show with an example how, still using spherical obstacles, the microscopic 
arrangement of the obstacles can alter the effective diffusion coefficient, that is, two 
structures with the same porosity may have different Dg. Specifically, we consider 
the hexagonal packing configuration in two dimensions. The periodic cell for such 
a configuration is w = [—1/2,1/2] x [—•\/3/2, •\/3/2] with one disk at the centre and 
a quarter of a disk at each corner (see Figure 2.3). We solve this new cell problem 
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£ = 0.2 £ = 0.4 £ = 0.5 



Fig. 2.3. Periodic cell uJv{'x.) of the hexagonal configuration in two dimensions for different 
values of the obstacles’ radius e. The third cell corresponds to close hexagonal packing. 


and evaluate the effective diffusion as a function of </> using (2.18b). The resulting 
effective diffusion tensor is again a multiple of the identity; we plot the scalar Dg as a 
function of solid volume fraction (f) in Figure 2.2(a). We find that, for the same value 
of there are differences in the effective diffusion coefficient in a square or hexagonal 
configuration, which become important for cj) above 50%. 

As mentioned in the introduction, the exact values of De((j)) in ordered lattices 
can be computed using Rayleigh’s multipole method, in the form of an infinite system 
of algebraic equations involving image multipoles and lattice sums [26]. By truncating 
the infinite system, Rayleigh obtained a closed-form approximation for a square lattice 
(S) and cubic lattice (C). The corresponding result for a two-dimensional hexagonal 
lattice (H) was derived in [18] as: 


d!{4>) 


1 


1 - 
1 

1 

IF 


1 - 


1 - 


1 - 


2(j) 


l + cj)- 0 . 3058^4 

2(j) 

1 + cf)- 0.0754206 

3(/) 

2 + (j)- 0.3914(()io/3 


(<^ < 0.7) 

(2.20a) 

{<i) < 0.8) 

(2.20b) 

((^ < 0.25) 

(2.20c) 


The values in parenthesis indicate the maximum volume fraction for which the closed 
forms are valid [18].^ These formulae are also plotted in Figure 2.2, and show good 
agreement with the multiple-scales result. 


2.2. Extension to non-periodic structures. Until now we have considered 
porous media with a simple regular structure, with the centers of the obstacles lying 
on a regular lattice, and accounted for macroscopic gradients in porosity by allowing 
the radius of the obstacles to vary slowly with position. In this section we consider 
an alternative way in which the porosity may vary: we suppose that the obstacles 
are all the same size, but that the number density of obstacles varies with position. 
This means that the obstacles can no longer be arranged on a regular periodic lattice. 
However, we suppose that their arrangement is locally a periodic lattice, but that the 
scale and orientation of the lattice vary slowly with position. To make this precise, we 
assume that the centers can be mapped to a regular periodic lattice with a map that 
depends only on the slow scale. Effectively this map defines (slow) curvilinear coor- 


^The exact result of the multipole method for various geometries in two and three dimensions 
is computed in [18] by direct inversion of the matrix, for a growing number of multipoles (with the 
matrix becoming larger and larger) until the result converges. 
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dinates in which the structure is periodic. We use a recently developed generalization 
of the multiple-scales method that can handle such a microstructure [27]. 

Thus we suppose there is exists a transformation W : 12 —^ 12' mapping the 
centers of the obstacles into a regular lattice with periodicity S, as illustrated in 
Figure 2.1. We could now apply the method of multiple scales in the transformed 
domain, averaging over the fast scale, before inverting the transformation to write the 
homogenized equation back in the original variables. However, since the fast and slow 
scales are treated as independent, there is no point in transforming the slow scale, 
only to invert the transformation later. 

Instead we suppose that the solution C is a function of the slow scale x and the 
transformed fast scale 


_ W(x) 

^ “ 6 ’ 

and that these variables are independent. Then, since the microstructure is periodic 
in y, we assume that C is periodic in y, with period one. We denote the unit cell in 
y by w and the obstacle by y G a;s(ri), where 

Ws(ri) = {y G w : (5y G W{B^{r,)) - W(ri)}. 


The available volume is then ujy (x) = w \ Wg (x). 

Using the chain rule spatial derivatives now transform according to 


d did 

dxi dxi S dyj ’ 


( 2 . 21 ) 


where Fij = dWjjdxi or F = and we are using the summation convention for 
repeated indices. 

The derivation of the homogenized equation is similar to the perfectly periodic 
case. The result is again an advection-diffusion equation for the concentration c as 
in (2.19), 


m 


= Vx- 


De{x)V^C 


£>e(x)VxV’(x) ■ 

V'(x) 


with a modified diffusion tensor 


De{x) = Id 


V’(x; 


-F(x) 



Jr (x, y) dy 


U-i(x). 


The function T = T^ now satisfies the following cell problem 


Vy • (F^FVyU) = 0 

F'^F VyTj • Hy = F'^F Hy ’ 6^ 

Fi periodic 


y G a;„(x), 
y G i9w„(x), 
in y. 


(2.22a) 


(2.22b) 


(2.23a) 

(2.23b) 

(2.23c) 


In general the cell problem (2.23) now depends on the slow scale x not only through 
cdy(x) but also because F depends on x. 
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2.2.1. Conformal maps. Since the Jacobian matrix F is dependent only on 
the slow scale x it represents a constant linear transformation on the fast scale. Such 
a transformation can always be written 

F = UDV 

where U and V are real orthonormal matrices and U is a diagonal matrix. In fact, 
in our application we expect U and V to have positive determinant and therefore 
be rotations (so that we maintain the right-handedness of the coordinate system). 
Thus F can be thought of as successive application of a rotation, a stretch along the 
coordinate axes, and then another rotation. 

In the particular case in which the stretch is isotropic, so that the entries along 
the diagonal in D are all equal, then F reduces to an isotropic stretch and a single 
rotation. In this case we may write F = aR, where a ^ 0 is a scalar stretch, and 
i? is a rotation matrix. Maps W(x) whose Jacobian F have this property preserve 
angles, that is they are conformal. Because the stretch is isotropic, spherical obstacles 
remain spherical under such a transformation, with the new radius of the original 
sphere centered at x being 


e(x) = i||W(x + e)-W(x)||. (2.24) 

For conformal maps the cell problem simplifies considerably, both because the 
sphere is mapped to a sphere, but also because 

= a^RR^ = a^Id. 

Thus the cell problem (2.23) reduces to our original cell problem (2.13), with the slow 
scale felt only through the change in obstacle radius with position via (2.24). 

In two dimensions, we may generate conformal maps by taking advantage of com¬ 
plex variables. If we write z = x + iy then any holomorphic function W : If C C —>■ C 
such that W'(z) yf 0 for z G If is a conformal map. An example that we will use in 
our simulations later is the function 

= W{z) = ^ (2.25) 

2i Z 

which is one-to-one and holomorphic except at z = 2. In particular, this defines a 
conformal map between U = [—1/2,1/2]^ and VF(f2), with inverse W~^{z) = 2— 1/z. 
In coordinate form, 


/ _ 2 - X / _ y 

"" ( 2 - a :)2 + y 2 y ( 2_^)2 + y 2 - 

In Figure 2.4 we plot the original domain 17 and the domain 17' = IF(I7). We see that 
in the transformed variables x' the centers of the obstacles lie on a regular square 
lattice, but the radii of the obstacles vary with position, that is, the structure is very 
similar to the one depicted in Figure 2.1. 

3. A random disordered porous medium. In this section we take a very 
different approach to model the same physical problem. We assume that the solid 
matrix is still formed of solid immobile non-overlapping spheres or disks of radius 
e at positions G 17, but now the are randomly distributed. The non-overlapping 
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(b) - 


y'o.o 


“■-0.5 0 0.5 

X 

Fig. 2.4. Original domain Q and the mapped domain Q' under the transformation (2.25). 
Parameters used are: 5 = 0.02, e = 0.01. 

constraint means that the Rj cannot be distributed independently, but we suppose 
that they are distributed identically, that is, the probability distribution function is 
invariant with respect to permutations of the particle labels 1 ,..., Ng. 

We suppose that there are mobile solute particles at positions Xi(t), i = 
1,... Njn, at time t, which undergo Brownian motion with constant diffusion coefficient 
Do, which we may set to unity as before through appropriate nondimensionalisation. 
In contrast to the previous section, we relax the assumption that the solute particles 
are point particles and allow for them to have a finite size. To keep things simple, we 
assume the solute particles are also spherical and of radius 1 and interact via 

hard-core elastic collisions with each other as well as with the obstacles. We will use 
the same approach as in our previous works [5, 6], for which we require that the total 
particle volume fraction is small so that three-particle interactions and higher form 
a negligible fraction of state space and can be ignored. In terms of the parameters 
introduced above, this means that Nmf-ti 1^1 = 1- 

Each solute particle evolves according to the overdamped Langevin stochastic 
differential equation 

dX, = y2dWi, (3.1) 

where the are Nm independent d-dimensional standard Brownian motions. We 
suppose that the initial positions X^ (0) are random and identically distributed (that 
is, the probability distribution is invariant with respect to a permutation of particle 
labels). 

3.1. The Fokker—Planck equation. Let P{x, r, t) be the joint probability den¬ 
sity of the Nm + Ng = N particles, where x = (xi, ..., xat^) and r — (Ri,..., RjvJ- 
We introduce the diV—dimensional position vector ^ = (x,r) with the coordinates 
of the mobile particles in the first dNm positions, and the obstacle positions in the 
remaining dNg positions. The probability density evolves according to the Fokker- 
Planck equation 

— = V|P, (f,r-)Gflf. (3.2a) 

Although this looks like the standard diffusion equation (I.la) we must remember that 
(1.1a) is solved in the low-dimensional physical space C O while equation (3.2a) is 


(a) 


y 0 


O o 


o 

’ o < 
O o o < 
° o o o o o o 

o o o o o o o 

O O o o O o o 

O O O od 

o O < 


o o 
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solved in the high-dimensional configuration space corresponding to N 

copies of the physical domain 12 minus the set of illegal configurations corresponding 
to particles overlapping, 

s.t. lie,-C,|| <£* + £,}, 

where = Cm for i < and = e otherwise. It is convenient to introduce e^m 
as the distance between centers at contact between an obstacle particle and a mobile 
particle, e„„, = e„, -I- e. 

On the collision surfaces 90^ we have the reflecting boundary condition 

0 = n-VsP, (3.2b) 

where n is the projection of the unit normal on the first dNm coordinates. 

Since P{x,r,t) is invariant with respect to permutations of the labels in r, the 
marginal density functions of P corresponding to fixing the position of one obstacle 
and integrating over the positions of the remaining obstacles are all identical, and 
given by 

s(x) = I P(x, r, 0)(5(ri — x) daidr. (3.3) 

This gives the probability of finding an obstacle in a given position, that is, it is the 
obstacle population density. Since |Pe| is the volume of one obstruction, the local 
volume concentration of obstacles is 


(f){x) = iVs|^e|s(x) = 4 'S(X), 


(3.4) 


where, as before, 4> is the average solid volume fraction. 

In order to compare the Fokker-Planck equation (3.2) to the model we have 
derived via multiple scales, we need an equation for the marginal solute density p(x, t) 
defined on the physical domain 17. As with the obstacles, since the Nm mobile particles 
are identically distributed, P is invariant to permutations of their labels and we can 
define the solute population density as 

p(x, t) = / P{x,r,t)S{x.i — x) dxdr. (3.5) 

inf 

Then the concentration of mobile particles is t) = A^mP(x, t) and the normalized 

concentration (used in §2) is c = Cm/ Jq Cmdx = p. 

The procedure we adopt to derive an equation for p from (3.2) is a systematic 
asymptotic expansion as Nge‘^ + Nmf-m 0. We use matched asymptotic expansions 
in configuration space, with an outer region in which particles are well-separated, 
and an inner region in which two particles are close together. In contrast with other 
approaches this systematic procedure does not require a closure assumption. 

The derivation is analogous to that presented in our previous work [5], and we 
omit the details. In fact, the problem studied here can be regarded as a particular 
case of the model for two species of diffusing and interacting particles presented in [5] 
if we formally set the diffusivity of one of the species to zero. Although the derivation 
should really be modified in this limit, the result is identical to equation (22a) in [5], 
which, in our notation, reads 


dp 

dt 


Vx ■ [l)e(x)VxP - v(x)p] 


(3.6a) 
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where the diffusion and drift coefficients are 

I)e(x) = 1 + {N^ - - iV«^ef^s(x), (3.6b) 

v(x) = ^ ^^^ ef^Vxs(x), (3.6c) 

for d = 2,3. The effective diffusion coefficient of solute particles has three compo¬ 
nents: the base molecular diffusion, an enhanced diffusion due to collective motion 
of finite size solute particles, and reduced diffusion due to excluded-volume interac¬ 
tions with the obstacle particles. The effective advection u(x) is due to the gradient 
in the distribution of obstacles; it indicates that the particles are advected towards 
regions with fewer obstacles. Note that there is no system-size expansion in (3.6): 
the equations are equally valid with small numbers of solute particles or obstacles 
(for example, the collective term disappears if there is only one solute particle, as it 
should). 

3.2. Model for point particles diffusing in a stochastic porous medium. 

The two approaches we have taken each have their strengths and weaknesses. The 
multiple-scales approach can handle arbitrary volume fraction of obstacles, but is 
limited to locally periodic structures and the diffusion of point particles of solute. 
The Fokker-Planck approach can handle arbitrary obstacle configurations and finite- 
size solute particles, but can only be reduced to a low-dimensional effective diffusion 
equation in the limit of small volume fraction. 

In order to compare the two approaches later, we set Cm = 0 in (3.6) to consider 
point solute particles. The effective diffusion and drift coefficients reduce to 


£)e(x) = 1 - 


$ 


-s(x) = 1 - 




(d-1)-- ^ {d-iy 

v(x) = -$Vxs(x) = -Vx(('(x), 


(3.7a) 

(3.7b) 


where $ = We plot Dg given by (3.7a) in Figure 2.2 with dashed 

lines. We observe that it agrees with the multiple-scales value Df. for small volume 
fractions (p (we will show this formally in §4.1). We note that, in the cubic obstacle 
configuration (d = 3), the asymptotic value is a good approximation to De{(j)) 
throughout the whole range of valid porosities. 

For point solute particles, the stationary density poo(x) takes a very simple form. 
Substituting (3.7) in (3.6a) and imposing no-flux boundary conditions, we find 


Poo(x) = K 





-i d-1 


(d-i) 


1 - (j){x) = V'(x), 


(3.8) 


where k is the normalization constant. This corresponds to the uniform measure in 
the available space, and is consistent with the stationary density found via multiple 
scales. 


3.3. Diffusion through moving obstacles is easier than through fixed 
obstacles. Before moving on to compare our different approaches, we briefly exam¬ 
ine how the effective transport properties of the solute particles change when the 
obstacles themselves are allowed to diffuse, with molecular diffusion Ds (in our di¬ 
mensionless setting this represents the ratio of the diffusion coefficient of the obstacles 
to that of the solute). This would be relevant in biological applications such as the 
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diffusion through biological tissues or the cytoplasm, where one is interested in the 
diffusion of small molecules through an environment containing large macromolecules 
[22]. These macromolecules constitute the “solid phase”, which is sometimes con¬ 
sidered as a porous structure since its dynamics are much slower than the smaller 
diffusing particles. Using the general model in [5] the coefficients in (3.7) change to 


£>e(x,t) = 1 - 


1 <(>(x,t) 

l + D,{d-l)^ 


v(x,t) = 


{d — 1) + dDg 
{d-l)(l + Ds) 


Vx(('(x,t), 


(3.9) 


for d = 2,3. Setting Hg = 0 in (3.9) recovers the expressions in (3.7) as expected. 
From (3.9) we see that the faster the obstacles move (the larger Dg), the less they 
impede the diffusion of the point solute particles, with their effect disappearing com¬ 
pletely as Dg —>■ oo. On the other hand, the larger Dg is, the larger the coefficient 
in front of the drift term v(x). This drift does not disappear in the limit, with 
V = —c?Vx</>/ {d — 1) as Dg —>• oo. Of course, since the obstacles are much larger than 
the solute particles, we would expect that they diffuse more slowly, i.e. that Dg < 1 
in any practical situation. 

4. Comparison between methods. In this section we compare the macro¬ 
scopic models for diffusion in a porous medium of variable porosity which we derived 
via multiple scales in §2 and using the Fokker-Planck approach in §3. As we men¬ 
tioned above, the Fokker-Planck approach can only be systematically reduced to a 
low-dimensional effective diffusion equation in the limit of small volume fraction. We 
observed in Figure 2.2 that the multiple-scales-derived diffusion coefficient seems nu¬ 
merically to asymptote to the Fokker-Planck-derived diffusion coefficient in this limit. 
In the next section we show that this is indeed the case, by considering the asymptotic 
solution to the multiple-scales model in the limit of low obstacle volume fraction $. 
We then compare our effective equations with each other and with direct numerical 
simulations in a variety of test problems. 

4.1. Model for an ordered porous medium with low porosity. We con¬ 
sider the model (2.18) in the limit of low volume fraction 4>. This means that the 
local volume fraction (/)(x) is also small (almost everywhere), corresponding to small 
relative obstacle radius e(x) 1. In this limit, the cell problem (2.13) can be solved 
explicitly. 

Since x, and hence e(x), are constants as far as the cell problem (2.13) is con¬ 
cerned, we can look for an asymptotic solution to (2.13) in terms of the small param¬ 
eter e(x). Consider, say, the first component of the vector function P, which satisfies 


Vyri = 0 yew„(x), (4.1a) 

Vyri-y = yi ||y||=e(x), (4.1b) 

Fi periodic in y. (4-lc) 

We use the method of matched asymptotic expansions, supposing that the unit cell 
LUyipc) can be divided into two regions: an inner region in which ||y|| ~ 0{e), and an 
outer region in which ||y|| ^ e. In the inner region, we set y = e(x)Y and define 
7i (x, Y) = Fi (x, y) to give 

Vy7i = 0 

VY 7 i-Y = eYi on ||Y|| = 1, 


(4.2a) 

(4.2b) 
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with a matching condition as ||Y|| —>■ oo. Expanding 7i(x, Y) = 7l°Hx,Y) + 
e:7f’(x , Y) + • • • gives that the leading-order inner solution 7^*^^ is simply a con¬ 
stant in Y, whence the leading-order outer solution is also constant. At first order in 
e we find 

V|.7f ^ = 0 (4.3a) 

Vvvf ^ • Y = Yi on II Y|| = 1. (4.3b) 


Using polar {d = 2) or spherical {d = 3) coordinates, we look for a solution to (4.3) 
of the form 7^^^ = f{x,R) cos9 {d = 2) or 7^^^ = f{x,R)sm9cos(p {d = 3), where 
Yi = Rcos9 when d = 2 and Yi = i?sin0cos(/3 when d = 3. We find that 

/(x,i?)=A(x)A+^H^^^, (4.4) 

for an unknown function A(x). Since the leading-order outer solution is constant in 
y, matching gives A(x) = 0. Thus 

T’(x.Y) = -i^^ + 7''(x). (4.5) 

Thus the first non-constant term in the outer expansion is 0{£‘^). Matching with the 
outer solution gives that 

El(x,y) - constant-H -, (4.6) 

(d- 1) ||y||'* 

as ||y|| 0. 

We can now use this asymptotic behavior to determine the outer solution at this 
order. However, it is possible to determine the effective diffusion coefficient, which 
is our primary aim, with the information we already have. Since the integrals we 
have to evaluate in (2.18b) are all derivatives with respect to some component of y, 
by integrating with respect to this component first we turn the volume integral over 
the unit cell ujy(x) into surface integrals over the exterior periodic boundaries and 
the interior boundary with the solid obstacle. The contributions from the exterior 
boundaries cancel due to periodicity, while on the interior boundary we can use the 
asymptotic solution (4.5). The result is 


/ 

J OJ 




T ^ 2TTe‘^{x) 

Jt dy = — 


Thus as e —>■ 0, De is a scalar multiple of the identity, equal to 

~' - ~ ‘ ' - i0Tr 


(4.7a) 


since = 1 ~ 4'i^) = 1 ~ 2(d — l)7re‘^/d ^ 1 at leading order. As expected, 

this result agrees with the asymptotic value (3.7a) obtained with the Fokker-Planck 
approach. 

The drift term in the multiple-scales homogenized equation (2.19) is v(x) = 
i7e(x)Vxl/'(x)/'!/;(x). Using (4.7a), we obtain 


Z)e(x)Vx((>(x) ^ d-1- (j) 

l-((i(x) (d-l)(l-(/)) 


Vx((>(x) ~ -Vx<('(x). 


(4.7b) 
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This asymptotic value also agrees with the drift obtained in the reduced Fokker- 
Planck model; see (3.7b). 

Finally we comment that the nature of the calculation above makes it clear that 
the configuration of the inclusions does not affect the diffusion coefficient at this order, 
since the dominant contribution comes from the solution of the cell problem in the 
inner region, which has no information about the position of the inclusion(s) in the 
unit cell. 

4.2. Numerical simulations. The aim of this section is to compare the two 
models for diffusion in porous media against each other as well as against numerical 
simulations for the full problem and stochastic simulations of the particle system. 

First, we consider the mean squared displacement of particles diffusing in two 
homogeneous porous media with the same porosity, namely a deterministic structure 
with a square lattice of obstacles and a random structure with obstacle configurations 
drawn from a uniform distribution with non-overlapping constraints. 

Second, we consider the spreading out of a localized initial concentration of par¬ 
ticles in porous media with gradients in porosity. Again we consider locally periodic 
structures accessible to the multiple-scales analysis and random structures with the 
same (ensemble) average porosity. 

All the simulations are made in a two-dimensional unit square domain with 
Nm point mobile particles and Ng hard-disk obstacles of constant radius e. When the 
porous structure is random with probability law s(x), a new configuration of obstacles 
is generated for every new run. 

4.2.1. Effective diffusion coefficient via the mean squared displacement. 

In this section we compare the diffusion coefficient computed from simulations of 
the discrete stochastic system to the effective diffusion coefficient obtained in the 
previous sections, either from multiple scales De (2.18b) if the obstacles are placed in 
a regular structure, or from the Fokker-Planck description De (3.7a) if the obstacles 
are randomly distributed with density s(x). In particular, we consider two porous 
media with uniform porosity $ (so that s(x) = 1 so that the drift is v = 0) and the 
same number Ng of obstacles: (i) a square lattice configuration, and (ii) a random 
uniform configuration of obstacles with non-overlapping constraints. 

The numerical value of the diffusion coefficient is obtained from the mean-square 
displacement, using the relation (r^(t)) = 2dDet as t —>■ oo. To evaluate the mean- 
square displacement, we run M = 1000 runs with Nm = 100 mobile particles, and 
compute (r2(t)) = ll^f^ ^ (0)lP- HereXf^(t) is the position 

of the ith particle in the k realization at time 

Our stochastic simulations are performed integrating Eq. (3.1) using an Euler- 
Maruyama scheme, with reflective boundary conditions between mobile particles and 
obstacles (cAl^) and periodic boundary conditions on the outer boundary (9S1). The 
reflective boundary conditions between the mobile point particles and the obstacles 
are implemented similar to as in [5], namely, the distance that a particle has trav¬ 
elled (illegally) inside an obstacle is reflected back into the domain To do that, 
we compute the point on the obstacle boundary where the particle penetrated, and 
compute a particle-wall elastic collision on that point. The integration timestep At 
must be chosen carefully so that virtually no collisions are missed. A convergence 


II Since every realization is done with Nm point particles, this is equivalent to averaging over 
particle trajectories and regenerating the solid matrix every Nm realizations). 
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study is shown in the Appendix and based on this, we have used At = 1.2434 • 10“® 
in the results presented below. 

We perform experiments with the periodic and random porous media at porosities 
$ = 10,20 and 30%. From Figure 2.2(a), we expect differences between the random 
and periodic porous media to become apparent from $ = 20%. The question is 
whether the discrepancy between Dg and De is real (due to the structure) or artificial 
(due to the nature of the asymptotic approximation in obtaining -Dg). 

We plot the mean-square displacement (r^(t)) against time in Figure 4.1(a) for 
the square lattice (solid lines with error bars) and the random structure (dashed lines 
with crosses and error bars). (The error bars indicate the 95% confidence interval, 
or ±1.96SD values, of each data point, and are barely discernible.) As expected, (r^) 
increases linearly with time, with a slope that decreases with $. 



(b) ■ 
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Fig. 4.1. (a) Mean squared displacement (r^) as a function of time t for diffusion of point 
particles in the presence of obstacles. The obstacles are arranged on a square lattice (error bars) 
or uniformly distributed (error bars with crosses) and at volume fraction <I> = 10,20,30%. The 
theoretical curves using Dg (Eq. (2.22h)) and Dg (Eq. (3.7a)) are shown with solid lines and 
dashed lines, respectively, (b) The same data replotted as (r^)/(4t), to demonstrate that by the 
final simulation time the numerical diffusion coefficient has converged. The data points inside the 
dashed black rectangle are used to compute the numerical value of diffusion. In both plots, the error 
bars indicate the 95% confidence interval. Parameters used for the simulations are: Nm = 100, 
M = 1000 (a new obstacle configuration was generated at every new run in the random case). 
At = 1.2434 ■ 10“®, e = 0.0126, and Ns = 200,400,600 for <I> = 0.1, 0.2, 0.3 respectively. 


The diffusion coefficient is given by limt_>.oo(i’^(t))/(4t). To check that we have 
run the simulation for long enough and rule out any anomalous transient diffusion, in 
Figure 4.1(b) we plot (r^(t))/(4t) [2]. This type of plot highlights any time dependence 
in the diffusion coefhcient [30]. We observe that, for all curves in Figure 4.1(b), the 
curves have slope 0 after the first t = 0.05 and hence have converged. To evaluate 
the diffusion coefficient, we average over the last At = 0.05 of each simulation (data 
points marked with a dashed rectangle in Figure 4.1(b)). From Figure 4.1(b) and the 
estimated values (data not shown), we note that: (i) For the square lattice case, theory 
(eq. (2.22b), solid lines) and simulation results (solid error bars) for Dg agree very 
well, as expected since Dg from (2.22b) is exact, (ii) The random media simulation 
results (error bars with crosses) agree well with the asymptotic value Dg (eq. (3.7a), 
dashed lines) for 4) = 0.1,0.2 but there is a significant discrepancy for $ = 0.3. (hi) 
The difference between the regular and random porous media is within the error bars 
for $ = 0.1 but becomes apparent for 4> = 0.2, 0.3. 

Interestingly, while the multiple-scales effective diffusion coefficient Dg does a 
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better job for the period medium (as we might expect), the Fokker-Planck effective 
diffusion coefficient D^. seems to do a slightly better job for the random medium. This 
was not at all obvious, since this coefficient is only the leading term in an expansion 
as $ —>■ 0 (while the multiple-scales coefficient is valid for all <f>). It seems that, for 
a given obstacle volume fraction $, random porous media may have a slightly lower 
diffusion coefficient than periodic ones. Because He is an asymptotic expansion, we 
can use the order of the next term to estimate the error. The next term in the 
asymptotic expansion of De is O ((2e)'^iV^, {2e)^Ns) [6], which gives 0.065 and 0.1459 
for $ = 0.2 and 0.3 respectively. The discrepancies in Figure 4.1(b) between the 
asymptotic values De and the simulation results for $ = 0.2 and 0.3 are 0.0103 
and 0.0257, respectively. In other words, the effective diffusion De does better than 
expected from the asymptotic error bounds. 

4.2.2. Diffusion in a gradient of porosity. For our second model comparison, 
we consider a porous medium with a non-uniform porosity. As before, we consider 
both a locally periodic array of obstacles of constant radius e, and a random array of 
obstacles giving the same local porosity. 

For the locally periodic structure we use the arrangement illustrated in Fig¬ 
ure 2.4(a). The obstacles have a fixed radius e = 0.01, and the average volume 
fraction of obstacle is $ = NsTtc^ = 0.059. To generate a random periodic medium 
with the same (ensemble) average local porosity we need to determine the probability 
density function of obstacle position s(x), which we have seen is related to the volume 
concentration of obstacles by 


cj){x) = <l)s(x). 


(4.8) 


This is easily found once we have determined the variable obstacle radius e(x) in the 
mapped multiple-scales domain (shown in Figure 2.4(b)). 

To determine e(x), consider one representative cell A(x) centered at x in the 
original domain fl. The cell A(x) is mapped to a square of side d centered at x' in 
n', which we denote A'(x'). The area of A(x) is then 


|A(x)| = 


dx = 


/A(x) 


/A'(x') 


det J^y^dx' = 


'A'(x') det Jw 


dx' = 


/A'(x') 


dx' 
11^014 ■ 


(4.9) 


Since the volume fraction in the cell is conserved, we have 


cj){K) 


|A(X)| 


7’'(fe(x))^ 

<52 


7r£r(x)^. 


(4.10) 


Thus knowing |A(x)| through (4.9) allows us to determine both (/)(x) (and therefore 
s(x)), and also 


e(x) 


e 

Trai’ 


(4.11) 


which we need in order to solve the multiple-scales cell problem (2.13). In Figure 4.2(a) 
we plot the density s(x) corresponding to the configuration shown in Figure 2.4(a). 
We see that the maximum local solid volume fraction is about 2.5 times the average 
solid volume fraction $. In Figure 4.2(a) we show one realization of obstacles ran¬ 
domly drawn from the corresponding probability distribution, with a non-overlapping 
constraint. 
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Fig. 4.2. (a) Density of obstacles s(x) corresponding to the configuration in Figure 2.4(a). (b) 
One realization of obstacles randomly drawn from the corresponding probability distribution, with a 
non-overlapping constraint. Parameters used are: 5 = 0.02, e = 0.01. 

We suppose that at t = 0 a drop of solute of radius e is placed centered at 
xq = (0.0392,0). We have chosen xg such that the drop does not intersect with any 
of the obstacles, that is, ||xo — r^H > 2e for * = 1,... ,Ns. Thus the initial condition 
is (^(x, 0) = l/(7re^) for ||x — xo|| < e, and zero otherwise. Of course, this initial 
condition does not satisfy the requirement of the multiple-scales method that it varies 
slowly with respect to the obstacle separation; nevertheless we expect that it will 
quickly spread into a function which does. 

In Figure 4.3 we illustrate the time evolution of the solute density for four dif¬ 
ferent models. In Figure 4.3(a) we show the (numerically calculated) true solution 
for the locally periodic distribution of obstacles. In Figure 4.3(c) we show the (nu¬ 
merically calculated) true solution for one realization of the random configuration of 
obstacles, distributed according to the density function s(x). In Figure 4.3(b) we 
show the solution of the effective equation (2.18) derived through the multiple-scales 
method, while in Figure 4.3(d) we show the solution of the (intrinsic version of the) 
effective equation (3.6) derived through the Fokker-Planck approach. In each case 
we show the intrinsic average c, since this allows direct comparison between the solu¬ 
tions of the effective equations and the solutions of the real problems.** We see that, 
perhaps counterintuitively, the maximum of the solute concentration initially moves 
to the right, towards the region of low porosity. This is because the localized source 
spreads out in all directions, but spreads out more in the high porosity region. The 
increased diffusion in the high-porosity region lowers the concentration more, giving 
the impression that the solute is moving towards the low-porosity region. 

At least visually, the two effective models capture the behavior of the true solu¬ 
tions, both for the regular structure and even for a single realization of the random 
structure. To make the comparison a little easier, we show in Figure 4.4 a slice 
along the line y = 0 for the true locally periodic solution, the multiple-scales-derived 
solution, and the Fokker-Planck-derived solution. We see that the agreement is re¬ 
markably good. 

While the intrinsic average c is the natural variable to compare with a particular 
realization of the microstructure, the volume average c is the natural variable to 
compare with an ensemble average over a random microstructure (and is also usually 
the natural variable in an application of the effective equations). While c converges 


**For the Fokker-Planck case, we recall that p = c given the normalization condition on c, and 
hence we identify the intrinsic average as c(x, t) = c(x, t)/^^ = p(x, t)/'il’, where '0(x) = 1 — ^(x). 
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Fig. 4.3. The concentration C(x, t) at times t = 0.1, 0.2, 0.3 for (a) a locally periodic distribution 
of obstacles; and (c) one realization of a random configuration of obstacles, distributed according to 
the same density function s(x). The intrinsic average c{x,t) at the same times from (b) equation 
(2.18) derived through the multiple-scales method; and (d) the equation (3.6) derived through the 
Fokker-Planck approach. The position of the initial solute drop is shown as a black disk. 



Fig. 4.4. The intrinsic average c(x, i) along the line y = 0 at times 1 = 0.1, 0.2,0.3, 0.4 computed 
from the multiple-scales equation (2.18) (blue, full) and from the Fokker-Planck equation (3.6) (red, 
dot-dashed). The true solution C for a locally periodic distribution of obstacles is shown in black 
(broken by obstacles). 


to the uniform measure as time evolves (as in Figure 4.3), the volume average c 
approaches a non-uniform density which is a multiple of the porosity. We show in 
Figure 4.5 the evolution of the volume average c for three models. In Figure 4.5(a) we 
show the solution of the effective equation (2.22) derived through the multiple-scales 
method, while in Figure 4.5(c) we show the solution of the effective equation (3.6) 
derived through the Fokker-Planck approach. In Figure 4.5(b) we show average of 
the true solution over 100 realizations of the randomly distributed obstacles. Again, 
visually, both models seem to do a good job of approximating the ensemble average. 
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Fig. 4.5. The volume average at times t = 0.1,0.2, 0.3 computed from (a) the effective 

equation (2.22) derived through the multiple-scales method; (h) the ensemble average of the true so¬ 
lution over 100 realizations of random configurations of obstacles distributed according to the density 
function s(x); and (c) the effective equation (3.6) derived through the Fokker-Planck approach. The 
position of the initial solute drop is shown as a black disk. 


In Figure 4.6 we show the solutions in Figure 4.5 along the strip y = [—3A/2, 3A/2] 
for A = 1/21. In Figure 4.6(a) we show the solutions of the multiple-scales-derived 
and Fokker-Planck-derived effective equations, along with the average over 100 real¬ 
izations of the true solution for a random distribution of obstacles (as in Figure 4.5). 
In Figure 4.6(b) we show some of the individual realizations along with the ensemble 
average, to give an idea of the variance. 




Fig. 4.6. The volume average c(x,t) along the strip y = [—3A/2,3A/2], computed by averaging 
the solution over bins of width A = 1/21 in x and 3A iny. (a) Solution at times t = 0.1, 0.2, 0.3, 0.4 
computed from the multiple-scales equation (2.18) (blue, full) and the Fokker-Planck equation (3.6) 
(red, dot-dashed), and an ensemble average over 100 realizations of the true solution for a random 
distribution of obstacles (black dots), (b) Ensemble average (black dots) at time t = 0.1 along with 
10 individual realizations of the true solution (colored), again averaged over bins of width A = 1/21 
in X and 3A in y. 
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5. Discussion. We have investigated the problem of diffusion through a porous 
medium in which the porosity is non-uniform. Transport through porous media is 
inherently a multiscale process, with properties determined by individual pores at 
the microscale while one is usually interested transport over much larger distances. 
To derive an effective macroscopic model, a common approach is to assume that all 
heterogeneities lie within the microscale and that, once these are averaged out, we are 
left with a homogeneous porous medium at the macroscale. It is tempting to suppose 
that, if the properties such as porosity do vary macroscopically, then all we have to do 
is carry out this procedure locally at each point (that is, treat the material as though 
it were uniform with a porosity equal to the local porosity), so that the effective 
diffusivity becomes a function of macroscopic position. However, in this paper we 
have seen that this is not the case, and that a more careful analysis is needed. 

We have considered two different approaches to the upscaling problem, suitable 
for deterministic and random porous media respectively, and have generalized them 
to heterogeneous media. The result is a macroscopic advection-diffusion equation, 
with the advection term accounting for the macroscopic gradients in porosity. 

First, we have extended the method of multiple scales to account for non-uniform 
porosity (via a microscopic cell geometry parametrized by the macroscopic variable) 
and non-periodic structures (providing they are locally periodic, that is, they can 
be mapped into periodic structures by a transformation depending only on the slow 
scale). The resulting equation is equivalent at leading order to the one presented in 
[32] using a volume-averaging approach, although some formally higher-order terms 
are included in their unit cell problem (which leads them to some worrying conclusions; 
for example, the effective diffusivity depends on the location of the centroid in the unit 
cell). The multiple-scales method has the advantage of being a systematic asymptotic 
expansion (for which higher-order terms could in principle be calculated) which is able 
to handle any porosity. Its disadvantage is that, even with our extensions, it requires 
the microstructure to be locally periodic. 

Our examples have considered spherical obstacles that remained spherical when 
mapped to a periodic arrangement, which has greatly simplified parts of the presen¬ 
tation. However, the technique is more general. In particular if the diffusion tensor 
is anisotropic (for example, if the inclusions were ellipsoids instead of spheres) the 
technique demonstrates systematically that the principal directions of the diffusion 
tensor would be aligned with the local axes of the ellipsoids. 

One question that always arises when deriving effective equations with the method 
of multiple scales is how much of an error is introduced by treating the microstructure 
as periodic, when in reality it is unstructured. To address this question, we considered 
a second approach involving diffusion through a random distribution of spherical 
obstacles. In this case we again used a systematic asymptotic expansion, but this 
time in the limit of low volume fraction of obstacles. The macroscopic equation turns 
out to be a particular case of our model for the diffusion of binary mixtures of finite¬ 
sized particles [5], when the diffusivity of one of the species is set to zero. 

We compared the macroscopic models described above to each other and also to 
their microscopic counterparts. For the examples we considered both effective material 
models performed well, and the differences between the structured and unstructured 
media were small. 

For both deterministic and random media, the mean-square displacement (r^) 
analysis in Figure 4.1 showed that the effective transport mechanism of the solute 
particle is still normal diffusion, that is, (r^) grows linearly in time, and it is only the 
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diffusion coefficient that changes from one (in the microscale) to in the macroscale. 
However, depending on the timescale used in such analysis, one can observe transient 
anomalous diffusion (the transient period can only just be seen in our simulations, see 
Figure 4.1(b)). In the context of porous media, if the ratio e between pore scale and 
typical domain length is not sufficiently small, diffusion may occur in a lengthscale 
shorter than the crossover region between transient anomalous diffusion and normal 
diffusion [30]. This transient behavior should not be confused with real stationary 
anomalous diffusion [2]. 

Appendix. Numerical convergence of mean-square displacement sim¬ 
ulations. Using the method described in §4.2.1, we compute the effective diffusion 
coefficient in porous media with 20% obstacle volume fraction, with obstacles either 
on a square lattice or uniformly distributed. We use the same parameters as in Fig¬ 
ure 4.1 except the timestep At, which we will vary to perform a convergence study. 
In addition, we use this study to check the number of trajectories (the number of par¬ 
ticles Nm times the number of runs M) required to obtain an accurate mean-square 
displacement. This is particularly important in the random case, as we need to obtain 
an accurate mean from the distribution of random media. 

The mean displacement of a diffusive particle according to (3.1) is h = \/2At. 
If the obstacle’s diameter is 2e, then we should use At such that h <g; 2e. We use a 
simulation time t = 0.25 (so that each particle has had time to diffuse across almost 
the whole domain fl). We choose a mean step h = (2e)/2^ for k = 0,1,... ,4, or, 
equivalently, a timestep At = where e = 0.0126. The diffusion coefficients 

are obtained as the average of the last 10 data points of the time series (after checking 
it has converged to a stationary value, see dashed area in Figure 4.1(b)). The resulting 
values for the periodic and random media, denoted by and respectively, are 
shown in Table A.l. The standard deviation in both cases is below 0.006 for all k. 

Table A.l 

Diffusion coefficient in a square lattice with volume fraction $ = 20% from simulations using 
At = Theoretical predictions De and De from Eqs. (2.22b) and (3.7a) respectively are 

also shown. 





Df'> 


Df'> 

De 

o 

II 

0.984483 

0.860079 

0.790490 

0.833896 

0.833425 

0.833163 

k = 1 

0.891180 

0.794840 

0.833218 

0.833427 



k = 2 

0.818925 

0.830819 

0.833424 




k = 3 

0.827845 

0.833261 





k = 4 

0.831907 












De 

o 

II 

0.895929 

0.815169 

0.779048 

0.811383 

0.807609 

0.8 

k = 1 

0.835359 

0.781305 

0.810877 

0.807624 



k = 2 

0.794819 

0.809029 

0.807675 




k = 3 

0.805476 

0.807759 





k = A 

0.810293 







The order of convergence of and D^e'^ as At —>• 0 is one (in At), as expected 
from the Euler-Maruyama integration scheme. We can apply repeated Richardson 
extrapolations and to improve the accuracy of these numerical results (see 
Table A.l). By the last extrapolation, appears to have five correct figures, three 
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^ (' 5 ') 

of which coincide with the theoretical value D^. Similarly, De has four correct fig¬ 
ures, two of which agree with the theoretical value D^. 
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